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We use a popular fictional disease, zombies, in order to introduce techniques used in modern 
epidemiology modelling, and ideas and techniques used in the numerical study of critical phenomena. 

We consider variants of zombie models, from fully connected continuous time dynamics to a full scale 
exact stochastic dynamic simulation of a zombie outbreak on the continental United States. Along 
the way, we offer a closed form analytical expression for the fully connected differential equation, 
and demonstrate that the single person per site two dimensional square lattice version of zombies 
lies in the percolation universality class. We end with a quantitative study of the full scale US 
outbreak, including the average susceptibility of different geographical regions. 
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I. INTRODUCTION 

Zombies captivate the imagination. The idea of a 
deadly disease that not only kills its hosts, but turns 
those hosts into deadly vectors for the disease is scary 
enough to fuel an entire genre of horror stories and films. 
But at its root, zombism is just that - a (fictional) disease 
- and so should be amenable to the same kind of anal¬ 
ysis and study that we use to combat more traditional 
diseases. 

Much scholarly attention has focused on more tradi¬ 
tional human diseases |1Q], but recently, academic atten¬ 
tion has turned a bit of thought onto zombies as a unique 
and interesting modification of classic disease models. 
One of the first academic accounts of zombies was the 
2009 article by Munz et al. 12], in which an early form 
of a compartmental model of zombism was introduced. 
Since then, there have been several interesting papers 
published including works that perform Bayesian esti¬ 
mations of the zombie disease parameters [22] . look at 
how emotional factors impact the spread of zombies m, 
using zombies to gain insight into models of politics [9], 
or into the interaction of a zombie epidemic and social 
dynamics mm- Additional essays can be found in two 
books collecting academic essays centered around zom¬ 
bism [41 [20] 

Besides the academic papers, zombies have seen a 
resurgence in fiction. Of particular note are the works 
of Max Brooks, including a detailed Zombie Survival 
Guide [I], as well as an oral history of the first zom¬ 
bie war ^ in a hypothesized post outbreak world. In 
both these works Brooks provides a rich source of infor¬ 
mation about zombies and their behavior. In particular, 
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he makes the connection to disease explicit, describing 
zombies as the result of a hypothetical virus, Solanum. 

Zombies form a wonderful model system to illus¬ 
trate modern epidemiological tools drawn from statistical 
mechanics, computational chemistry, and mathematical 
modeling. They also form an ideal vehicle for public out¬ 
reach: the Center for Disease Control uses preparation 
for a zombie apocalypse mm to promote emergency 
preparedness. In this work, we will build up to a full- 
scale simulation of a zombie outbreak in the continen¬ 
tal United States, with realistic values drawn from the 
literature and popular culture (section [V| simulation ac¬ 
cessible online HU)- Before that, we shall use statistical 
mechanics to scrutinize the threshold of zombie viru lenc e 
that determines whether humanity survives (section |lV| ). 
Preceding that, we shall show how methods from com¬ 
putational chemistry can be used to simulate every indi¬ 
vidual h eroi c encounter between a human and a zombie 
(section [TTT|). But we begin by describing and analyzing 
a simple model of zombies (the SZR model) - the sim¬ 
plest and most natural generalization to the classic SIR 
(Susceptible-Infected-Recovered) model used to describe 
infectious disease spread in epidemiology. 


II. SZR MODEL 

We start with a simple model of zombies, the SZR 
model. There are three compartments in the model: S 
represents the susceptible population, the uninfected hu¬ 
mans; Z represents the infected state, zombies; and R 
represents our removed state, in this case zombies that 
have been terminated by humans (canonically by destroy¬ 
ing their brain so as to render them inoperable). There 
are two transitions possible: a human can become in¬ 
fected if they are bitten by a zombie, and a zombie can 
be destroyed by direct action by a human. There are 
two parameters governing these transitions: /3, the bite 
parameter determines the probability by which a zombie 
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will bite a human if they are in contact, and n the kill 
parameter that gives the probability that a human kills 
the zombie. Rendered as a system of coupled differential 
equations, we obtain, for a particular interaction site: 


Given the analytical solution, it is clear to see that the 
sign of P governs whether there will eventually be hu¬ 
mans or zombies in the final state. If a < 1,P > 0, 
so 
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Notice that these interactions are density dependent , in 
the sense that the rate at which we convert humans to 
zombies and kill zombies is dependent on the total count 
of zombies and humans in this site. This is in contrast 
with most models of human diseases, which frequently 
adopt frequency dependent interactions wherein S, Z, R 
would have been interpreted as the fraction of the popu¬ 
lation in the corresponding state. 

This distinction will become stark once we consider 
large simulations with very inhomogeneous populations. 
By claiming that zombies can be modeled by a single bite 
parameter /? that itself is a rate per person per unit time, 
we are claiming that a zombie in a block with 5,000 peo¬ 
ple would be one hundred times as effective at infecting 
new zombies as a zombie in a block with fifty people; 
similarly the zombie in question would be killed one hun¬ 
dred times faster. This would seem false for an ordinary 
disease like the flu, but in the case of zombies, we argue 
that it is appropriate. Zombies directly seek out hosts to 
infect, at which point the human and zombie engage in 
a duel to the (un)death. 

To facilitate analysis we can nondimensionalize the 
equations by choosing a relevant population size TV, and 
recasting in terms of the dimensionless time parameter 
r = t/3N and dimensionless virulence a = k//3 
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Unlike a traditional disease (e.g., as modeled by SIR), 
for the zombie model, we have a stable configuration 
when either the human or the zombie population is de¬ 
feated (S = 0 or Z = 0). Furthermore, unlike SIR , SZR 
admits an analytical solution, assuming R( 0) = 0, and 
with Zq = Z(0), So = S'(O): 


and the system will always flow to a final state composed 
of entirely zombies and no humans, where P denotes the 
number of zombies that survive. 

If however, a > 1, humans are more effective at killing 
zombies than zombies are at biting humans. With enough 
zombies in the initial state, we can still convert all of the 
humans before they have time to kill all of the zombies. 

We can recast the dynamics in terms of the variables 
P = Z P {1 — a)S and y = S/Z to gain further insights. 
First note that: 
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Hence if we choose TV = |P|, we end up with the very 
simple dynamics: 


P\r)= 0 (19) 
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( 20 ) 
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(5) Here we see that the dynamics is simply an exponential 
decay or increase in the ratio of humans to zombies y = 

(6) S/Z. The final populations in either case are easy to see 
due to the conservation of P. If zombies win we have 

(7) Z oo = Z 0 + (l- a)S 0 (24) 

(8) And if humans win 

(9) Soo = So - (25) 

a — 1 
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1. SIR model 

This dynamics should be compared to the similarly 
nondimensionlized density-dependent SIR model: 
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cause in the SIR model our infected population recovers 
on its own. This is contrasted with SZR , where the pro¬ 
cess of infection and recovery have the same functional 
form, depending on the product SZ. This /i is the inverse 
of the usual Rq parameter used to denote the infectivity 
of the SIR model, here used to make a closer analogy 
to the SZR model. It is this parameter that principally 
governs whether we have an outbreak or not. Unlike the 
a parameter for SZR which depends only on our dis¬ 
ease constants /3, /c, the relevant virulence for the density 
dependent SIR model (/i) has a population dependence. 

Notice again that while the only stable configuration 
for the SIR model is when there is no infected population 
(/ = 0), the SZR model is stable when either the humans 
or zombies are depleted (S = 0 or Z = 0). 

The SIR model does not admit a closed form analyt¬ 
ical solution, but we can find a parametric solution by 
dividing the first equation by the third, revealing. 

(H(t)-Rq) 

S(r) = Soe-(29) 

Using the observation that in the limit of infinite time, 
no infected population can persist, we can choose N to 
be the total population 

5o + /o + Ro = N = Sqq + i?oc (30) 

and so obtain a transcendental equation for the recovered 
population at long times. 

( Roo-Ro ) 

Roo=N- S 0 e (31) 

Unlike the SZR model, here we see that no matter how 
virulent the disease is, the epidemic will be self-limiting, 
and there will always have some susceptibles left at the 
end of the outbreak. This is a sharp qualitative differ¬ 
ence between zombies and more traditional SIR models, 
arising from the fact that the “recovery” of zombies is 
itself dependent on the presence of susceptibles. 

To visually compare the difference, in Figure[l]we have 
shown deterministic trajectories for both SIR and SZR 
for selected parameter values. 

III. STOCHASTIC SIMULATION 

While most previous studies modeling zombie popula¬ 
tion dynamics have been deterministic, things get more 



FIG. 1. Deterministic trajectories for the SIR and SZR mod¬ 
els with an initial population of 200 people, 199 uninfected 
and 1 infected. The (susceptible, infected, removed) popula¬ 
tion is shown in (blue, red, black) (color online). The SZR 
results are solid lines while the SIR results are lighter lines. 
For both models r = t/3N where N was taken to be the to¬ 
tal population. For the SZR model a was chosen to be 0.6, 
while for the SIR model /i was chosen to be 0.6 to show sim¬ 
ilar dynamics. Notice that in this case, in SZR the human 
population disappears and only zombies remain in the end, 
while the SIR model is self-limiting, and only a fraction of 
the population ever becomes infected. 



FIG. 2. Example Gillespie dynamics for the SIR and SZR 
models with the same parameter settings as Figure [l] The 
(susceptible, infected, removed) population is shown in (blue, 
red, black) (color online). The SZR results are solid lines 
while the SIR results are lighter lines. The two simulations 
were run with the same seed so as to match their dynamics 
at early times. 


interesting when we try to model discrete populations. 
By treating the number of zombies and humans as contin¬ 
uous variables in the last section, we are ignoring the ran¬ 
dom fluctuations that arise in small populations: even a 
ferociously virulent zombie infestation might fortuitously 
be killed early on by happy accident. Similar problems 
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arise in chemical reactions: reactions involving two types 
of proteins in a cell can be described by chemical reac¬ 
tion kinetics evolving their concentrations (like our SZR 
equations [4| , but if the number of such proteins is small, 
accurate predictions must simulate the individual binary 
reactions (each zombie battling each human). Interpret¬ 
ing our SZR transitions as reaction rates, gives us a sys¬ 
tem akin to a chemical reaction with two possible tran¬ 
sitions: 


(S,Z)^>(Z,Z) 
( S , Z) —► (S, R) 
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in final state 


When a human and zombie are in contact, the prob¬ 
ability of a bite in a small period of time is given by 
the bite rate and the size of the populations of the two 
species (/ 3SZdt ), and similarly for the probability of a 
kill. In order to efficiently simulate this dynamics, we 
use the Gillespie algorithm [7], which efficiently uses the 
computer to sequentially calculate the result of each one- 
on-one battle. 

The stochasticity gives more character to the simula¬ 
tion. The fully connected continuous dynamics modeled 
by the differential equation is straightforward: either the 
humans win and kill all of the zombies, or the zombies 
win and bite all of the humans. While the continuous ap¬ 
proximation may be appropriate at intermediate stages 
of the infection where the total population is large and 
there are a non-trivial number of infected individuals, 
we will eventually be interested in simulating an actual 
outbreak on an inhomogeneous population lattice, where 
every new site will start with a single infected individ¬ 
ual. But even though we may be interested in modeling 
the outbreak case (a < 1), we would like to allow the 
possibility that the humans manage to defeat the out¬ 
break before it really takes off. The stochastic Gillespie 
dynamics allows for this possibility. 

In Figure [2] we have shown an example of a single 
stochastic simulation using the same parameter settings 
as those used in Figure^ The stochastic trajectory over¬ 
all tracks the analytic result, but at points in the simula¬ 
tion there may be more or fewer zombies than anticipated 
if the dice fall that way. 

Another implication of stochastic dynamics is that it 
is not always guaranteed that a supercritical (a < 1) 
outbreak will take over the entire susceptible population. 
For the parameter settings used in Figure [l] and [2j namely 
a = 0.6 with a population of 200 and one infected indi¬ 
vidual to start, the zombies win only 40% of the time. 
Additionally, the number of zombies we end with is not 
fixed, as shown in Figure [3j 

In fact, we can solve exactly for the probability ^ext 
that an a < 1 simulation will go extinct in the limit 
of large populations, using an argument drawn from the 
traditional SIR literature. At the very beginning of the 
simulation, there is only one zombie, who will be killed 
with probability k/(/3 + k). If the first zombie is killed 


FIG. 3. Distribution for final zombies over 100,000 stochastic 
trajectories with the same parameters as Figure [2] Not pic¬ 
tured are the 60% of runs that end with no zombies in the final 
state. Compare these to the analytical result, in which the 
final population of zombies would be 81 with no possibility of 
surviving humans. 



FIG. 4. The observed fraction of simulations that end in 
an extinction for the zombie outbreak, f or 1 ,000 runs of 10 4 
individuals at various values of a (eqn. [33] ) . The observed 
extinction probabilities agree with the expectation that they 
should go as a , here shown as the dashed line. This is the 
same behavior as the SIR model. 


before it bites anyone, we guarantee extinction. Other¬ 
wise, the zombie will bite another human, at which point 
there will be two independent zombie lines that need to 
be extinguished, which will occur with probability P e 2 xt . 
This allows us to solve: 


PpTK t - 


AT 


P 


P + K 


p 2 
■*ext 


Text, — ^ — OL . 


(32) 

(33) 


The probability of extinction is just given by our dimen¬ 
sionless inverse virulence a. In Figure [4] we have shown 
the observed extinction probabilities for 1,000 Gillespie 
runs of a population of 10 4 individuals at various values 
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of a , and overlaid our expected dependence of a. 

This same extinction probability (P ex t = /jl = Rq 1 ) is 
observed for the SIR model m- This is not a coinci¬ 
dence. In precisely the limit that is important for study¬ 
ing the probability of an extinction event, namely at early 
times with very large populations, the SZR model and 
SIR are effectively the same, since the population of sus- 
ceptibles (S) is nearly constant. Writing S as So — SS, 
we have: 


dZ 
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Here as SS 0, the two models are the same with a = 
fiN/So, another indication that the density dependent 
SIR model’s virulence is dependent on population size. 

To get a better sense of the effect of the stochasticity, 
we can look at the mean fractional population in each 
state for various settings of a and choices for initial pop¬ 
ulation size. The results are shown in Figure [5] 
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FIG. 5. Mean final states as a function of model param¬ 
eters. One thousand different simulations are run for each 
cell. Each simulation starts with a single zombie or infected 
individual. The runs are run until they naturally terminate, 
either because the susceptible population is deleted, the zom¬ 
bie population is gone, or there are no more infected individ¬ 
uals. Each cell is colored according to the mean fraction of 
the population occurring in each state. The top row is for 
SZR simulations and the bottom row is for SIR simulations. 
In both cases N is chosen to be 100. Here the sharp contrast 
between density-dependent SZR and SIR is made apparent. 
Notice that density-dependent SIR is very strongly popula¬ 
tion dependent. 

Plotted are the fractional populations in the final state 


left for both the SZR model (top row) and SIR model 
(bottom row) for different parameter combinations of a 
and the initial population. In all cases, the N parameter 
was chosen to be 100. For each box, 1,000 independently 
seeded stochastic trajectories were calculated until com¬ 
pletion. Looking at the SZR results in the top row, we 
can see that the dynamics is fairly independent of popu¬ 
lation size once the population size gets above around 100 
individuals. The population dependence for lower popu¬ 
lation sizes is an effect of the stochasticity. We can clearly 
see a transition in the susceptible population near a = 1 
corresponding to where our continuous dynamics would 
show a sharp boundary. Here the boundary is blurred, 
again due to the stochasticity. The final dead zombie 
population R remains small for all values of a; for ex¬ 
tremely virulent zombies a 1, very few will be killed 
by the humans before all of the humans are converted, 
while in the other extreme few zombies are created so 
there are few to be killed. 

Contrast these results with the density dependent SIR 
dynamics shown in the second row. There can be no in¬ 
fected individuals left in the end, so only the fraction of S 
and R in the final state are shown. The two transitions in 
SIR couple differently to the population of infected and 
susceptible. While our nondimensionalized SZR model 
has Z' = (1 — a)SZ/N , our nondimensionlized SIR has 
I r = (S/N — fjb)I. This creates a very strong population 
dependence. The transition observed in the S population 
is largely independent of /i, except on the very small end. 
When we move to inhomogeneous population lattices this 
means that for the density dependent SIR model, the 
most important parameter governing whether a particu¬ 
lar site has a break-out infection is the population of that 
site on the lattice. 


IV. CRITICAL BEHAVIOR OF LATTICE 
MODEL 

Until now, we have considered fully connected, well- 
mixed populations, where any infected individual can 
infect any susceptible individual with equal probability. 
But surely, a zombie in New York cannot bite someone 
in Los Angeles. Investigation of the spatial spread of in¬ 
fectious diseases is an important application of network 
science ; social diseases spread among intimate contacts, 
Ebola spreads by personal contact in a network of care¬ 
givers, influenza can be spread by direct contact, through 
the air or by hand-to-mouth, hand-to-eye or hand-to-nose 
contact after exposure to a contaminated surface. For 
most diseases, dong bonds’ dominate the propagation to 
distant sites E2; airplane flights take Ebola to new con¬ 
tinents. Zombies do not fly airplanes, so our model is 
closer in spirit to the spread of certain agricultural infes¬ 
tations, where the disease spreads across a lattice of sites 
along the two-dimensional surface of the Earth (although 
not in those cases where pathogens are transported long 
distances by atmospheric currents). 
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To begin, we will consider a two-dimensional square 
lattice, where each site contains a single individual. Each 
individual is allowed to be in one of three states: S', Z, or 
R. The infection spreads through nearest neighbor bonds 
only. That is, a zombie can bite or be killed by any 
susceptible individuals in each of the four neighboring 
sites. 

To make direct contact with our zombie model, the 
rate at which an susceptible cell is bitten is given by /3Z 
where Z is the number of zombie neighbors (since S is 
one), and the rate at which a zombie site is killed is kS 
where S is the number of susceptible neighbors. 

Because all state transitions in the SZR model de¬ 
pend only on Z-S contacts, for computational efficiency, 
we need only maintain a queue of all Z-S bonds, that is 
connections along which a human and zombie can in¬ 
teract. At each step of the simulation, one of these 
Z-S bonds is chosen at random, and with probability 
/3/(/3 + k) = 1/(1 + a), the human is bitten, marking it 
as a zombie. We can then query its neighbors, and for 
all of them that are human, we can add a Z-S link to 
our queue. With probability k/(/3 + k) = a/(l + a) the 
zombie is killed, removing any of its links to neighbor¬ 
ing humans from the queue. This process matches the 
stochastic dynamics of our zombie model operating on 
the lattice. 

Simulating zombie outbreaks on fixed lattices, there is 
qualitatively different behavior for small a and large a. 
When a is large, the zombies do not spread very far, al¬ 
ways being defeated by their neighboring humans. When 
a is very small, the zombies seem to grow until they in¬ 
fect the entire lattice. This suggests evidence of a phase 
transition. Technically, the presence of a phase transi¬ 
tion would mean that if we could simulate our model on 
an infinite lattice, there should be some critical a (a c ), 
above which any outbreak will necessarily terminate. Be¬ 
low the critical value, there is the possibility (assuming 
the infection does not die out) of having the infection 
grow without bound, infecting a finite fraction of indi¬ 
viduals in the limit that the lattice size becomes infinite. 
The SIR model has been demonstrated to undergo such 
a phase transition, and we expect the zombie model does 
as well. 

The study of critical phenomena includes a series of 
techniques and analyses that enable us to study the prop¬ 
erties of phase transitions even on finite lattices. A major 
theme of critical phase transitions is the importance of 
critical points - where a system is tuned (here by varying 
a) to a value separating qualitatively different behaviors 
(here separating low-infectivity transient zombie infesta¬ 
tions from a potentially world-spanning epidemic). At 
critical points, the system can show scale free behavior ; 
there is no natural length scale to the dynamics, and 
various physical parameters will usually be governed by 
power laws (see below). 

With a chosen to be precisely at the critical value, 
we indeed see a giant component with fractal structure 
(Fig# Note that there are holes (surviving pockets of 


humans) of all sizes in the figure. This reflects the prox¬ 
imity to the threshold: the battle between zombies and 
humans is so evenly matched, that one gets an emergent 
scale invariance in the survival patterns. This is in keep¬ 
ing with studies of the SIR model, which shows a similar 
critical behavior and phase transition [Hj. 



FIG. 6. Example cluster resulting from the single population 
per site square lattice zombie model with periodic boundary 
conditions near the critical point a c = 0.437344654(21) on a 
lattice of size 2048 x 2048. 

Systems near critical points with this kind of scale 
invariance fall into universality classes. Different sys¬ 
tems (say, a real disease outbreak and a simple com¬ 
putational model) can in many ways act precisely the 
same on large scales near their transitions (allowing us to 
predict behavior without knowing the details of zombie- 
human (anti)social interactions). The SIR model on a 
two-dimensional lattice with a single person per site falls 
into the percolation universality class [3], though details 
of its cluster growth can differ [21]. Given that the SZR 
model has two second order couplings, it is of interest 
whether it falls into the same percolation universality 
class. 

To extract the scaling behavior of our zombie infesta¬ 
tion, we study the distribution P(s,a), the probability 
that a single zombie will generate an outbreak of size s 
at inverse virulence a. (An outbreak will be a fractal 
cluster in two dimensions, with ragged boundaries if it 
dies out before reaching the entire world.) At a = a c 
where the zombies and humans are equally matched, we 
have an emergent scale invariance. A large outbreak will 
appear to almost stop several times - it can be viewed 
as a sequence of medium-sized outbreaks triggering one 
another. Medium-sized outbreaks are composed of small 
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outbreaks, which are in turn composed of tiny outbreaks. 
At threshold, each of these scales (large, medium, small) 
is related to the lower scale (medium, small, tiny) in the 
same fashion. Let us oversimplify to say that at critical¬ 
ity an outbreak of size 3s is formed by what would have 
been three smaller outbreaks of size s which happened 
to trigger one another, and these in turn are formed by 
what would have been three outbreaks of size s/3. If the 
probabilities and form of this mutual triggering is the 
same at each scale, then it would not surprise us that 
many properties of the outbreaks would be the same, af¬ 
ter rescaling the sizes by a factor of three. In particular, 
we expect at the critical point to find the probabilities of 
outbreaks of size s to be related to the probabilities at 
size 5/3 by some factor /: 

P(S, «c) = fP(s/ 3, a c ). (36) 

This formula implies that P(s,a c ) oc s -r , with r = 
log(l//)/log(3). The distribution of epidemic infection 
rates is a power law. 

Figure [7] shows a thorough test of this dependence for 
our zombie model, following a procedure akin to that 
of reference m- We simulated a zombie outbreak on 
a two-dimensional lattice with periodic boundary condi¬ 
tions starting with a single zombie. With the outbreak 
sizes following a power law distribution, the probability 
that a site belongs to a cluster of size n s is P s = sn s , so 
that at the critical point P s ~ s 1_r . Integrating from s 
to 00 , the probability that a point belongs to a cluster 
of at least s in size (P> s ) should at the critical point it¬ 
self follow a powerlaw: P> s ^ s 2 ~ r . To find our critical 
point <a c , we ran many simulations until our integrated 
cluster size distribution followed a power law, using the 
interpolation methods of reference [21 to get a precise 
estimate of the critical point. 

For zombies on a two dimensional lattice, this critical 
point occurs at a c = 0.437344654(21), the resulting in¬ 
tegrated cluster size distribution is shown at the top of 
Fig. 0 Percolation theory predicts r = 187/91 in two di¬ 
mensions, and we test that prediction in the bottom part 
of Fig. [7| Here, if we were precisely at the critical point 
and the SZR model is in the percolation universality 
class, with infinite statistics we would have asymptoti¬ 
cally a perfectly straight line. Notice the small vertical 
scale: our fractional fluctuations are less than 0.1%, while 
our experimental results vary over several order of mag¬ 
nitude. The clear agreement convincingly shows that the 
zombie model on the two dimensional lattice is in the 
percolation university class. 

As an additional check, we computed the fractal di¬ 
mension of our clusters near the critical point using box 
counting, a distribution for which is shown in Figure [8] 
We find a fractal dimension D = 1.8946(14), compared 
to the exact percolation value of D = 91/48 = 1.895833. 

Why did we need such an exhaustive test (many 
decades of scaling, many digits in our estimate of a c )? 
On the one hand, a much smaller simulation could have 
told us that there was emergent scale invariance and frac- 




FIG. 7. The cumulative distribution of epidemic sizes for 
the two dimensional zombie model near the critical virulence. 
The critical point found was a c = 0.437344654(21). The top 
plot shows the probability of a site being in a cluster of at 
least s in size (P> s ). The fact that it forms a straight line 
on a log-log plot indicates that P> s is a power law, and the 
slope is 2 — r. For comparison, the red (color online) line 
shows the powerlaw corresponding to the percolation critical 
exponent: r = 187/91. The bottom plot shows the same data 
times s T ~ 2 using the exponent from percolation theory. The 
plot is very nearly flat suggesting the percolation exponent 
accurately describes the zombie model. 


tal behavior near the transition; one or two decades of 
scaling should be convincing. But it turns out that there 
are multiple different universality classes for this kind of 
invasion process, and their exponents r and D are rather 
similar. And a small error in a c can produce large shifts 
in the resulting fits for r and D - demanding efficient 
programming and fast computers to achieve a definitive 
answer. 

We conclude that the single person per site zombie 
infestation, near the critical virulence, will on long length 
scales develop spatial infestation patterns that are well 
described by two-dimensional percolation theory. 
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FIG. 8. A histogram of the observed fractal dimension of the 
zombie epidemic clusters as measured by box counting. These 
give a measured value of D — 1.8946(14), consist with the 
exact percolation fractal dimension of D — 91/48 = 1.895833. 


V. US SCALE SIMULATION OF ZOMBIE 
OUTBREAK 

Having explored the general behavior of the zombie 
model analytically, stochastically and on homogeneous 
single person lattices, we are prepared to simulate a full 
scale zombie outbreak. 


A. Inhomogeneous Population Lattice 

We will attempt to simulate a zombie outbreak occur¬ 
ring in the United States. This will be similar to our lat¬ 
tice simulation, but with an inhomogeneous population 
lattice. We based our lattice on code available for creat¬ 
ing a “dot map” based off the 2010 US Census data Q51- 
The 2010 Census released census block level data, detail¬ 
ing the location and population of 11,155,486 different 
blocks in the United States. To cast these blocks down 
to a square grid, we assigned each of the 306,675,005 re¬ 
ported individuals a random location inside their corre¬ 
sponding census block, then gridded the population into 
a 1500 x 900 grid based on latitude and longitude co¬ 
ordinates. The resulting population lattice can be seen 
in the top half of Figure [9] You will see the presence 
of many empty grids, especially throughout the western 
United States. This disconnects the east and west coasts 
in a clearly artificial pattern - our zombies in practice 
will gradually wander through the empty grid points. To 
add in lattice connectivity, we did six iterations of binary 
closing (an image processing technique) on the popula¬ 
tion lattice and added it to the original. The effect was 
to add a single person to many vacant sites, taking our 
total population up to 307,407,336. The resulting popu¬ 
lation map is shown in the bottom half of Figure [9] This 
grid size corresponds to roughly 3 km square boxes. The 
most populated grid site is downtown New York City, 
with 299,616 individuals. The mean population of the 


occupied grid sites is 420, the median population of an 
occupied site is 13. 


io 5 



FIG. 9. A 1500 x 900 grid of the 2010 US Census Data. The 
above figure gives the raw results. Notice the multitude of 
squares with no people in them in the Western United States. 
The bottom figure shows the resulting map after 6 steps of 
binary closing added to the original population. 


B. Augmented Model 

In order to more ‘realistically’ simulate a zombie out¬ 
break, we made two additions to our simplified SZR 
model. The first was to add a latent state E (Exposed). 
The second was to introduce motion for the zombies. 
Considered as a system of differential equations, we now 
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have: 


Si = -/3SiZi 

(37) 

Ei = —uEi 

(38) 

Zi = nEi uSiZi 

(39) 

R% — KjSiZi 

(40) 

II 

M 

1 

(41) 

<i> 


set of reactions: 


(Si, ZEE (£. -l,Ei + 1) 

(42) 

.. rp 

C Zi,Ei ) - (Zi + l,Ei-l) 

(43) 

k S-7- 

(Zi,Ri) * *> (Zi-l,Ri + l) 

(44) 

(ij) : (Z h Zj) ( Z . _ i, Zj + 1) . 

(45) 


Here i denotes a particular site on our lattice, (j) denotes 
a sum over nearest neighbor sites, (ij) denotes that i and 
j are nearest neighbors. In this model, zombies and hu¬ 
mans only interact if they are at the same site, but the 
zombies diffuse on the lattice, being allowed to move to 
a neighboring site with probability proportional to their 
population and some diffusion constant (/i). We assume 
that the humans do not move, not only for computa¬ 
tional efficiency, but because, as we will see, the zombie 
outbreaks tend to happen rather quickly, and we expect 
large transportation networks to shut down in the first 
days, pinning most people to their homes. The addition 
of a latent state coincides with the common depiction 
that once a human has been bitten, it typically takes 
some amount of time before they die and reanimate as a 
zombie. If a human is bitten, they transition to the E 
state, where at some constant rate (y) they convert into 
the zombie state. 

To choose our parameters we tried to reflect common 
depictions of zombies in movies. The work of Witkowski 
and Blais [22] performed a Bayesian fit of a very similar 
SZR model to two films, Night of the Living Dead , and 
Shawn of the Dead. In both cases, the observed a was 
very close to 0.8. This means that the zombies in the 
films are 1.25 times more effective at biting humans than 
the humans are at killing the zombies. We will adopt this 
value for our simulation. For our latent state, we adopt 
a value close to that reported for Shawn of the Dead , 
namely a half-life of 30 minutes. To set our movement 
parameter, we estimate that zombies move at around 1 
ft/sec. To estimate the rate at which the zombies will 
transition from one cell to the next, we assume that the 
zombies behave like a random gas inside the cell, so that 
the probability that a zombie will cross a cell boundary 
is roughly ^-^LvAt, that is, one-fourth of the zombies 
within vAt of the edge will move across that edge in 
a small amount of time. This suggests a value of /a of 
0.0914 /hr. This corresponds to an average time between 


P 

3.6 x 10 3 /hr/person 

a 

0.8 

n 

a/3 

T) 

2 /hr 

L 

0.0914 /hr 


TABLE I. The parameters chosen for our US-scale simula¬ 
tions of a zombie outbreak. These parameters were chosen to 
correspond with standard depictions of zombies and simple 
physical estimations explained in the main text. 

transitions of around 11 hours, which for a zombie stum¬ 
bling around a 3 km block agrees with our intuitions. 
Finally, to set a rate for our bite parameter, we similarly 
assume that the zombies are undergoing random motion 
inside the cell at 1 ft/sec, and they interact with a hu¬ 
man anytime they come within 100 feet. We can then 
estimate the rate at which humans and zombies will in¬ 
teract as SZ ^fjt 1 , which corresponds to a choice of (3 of 
around 3.6 x 10“ 3 /hr. Another way to make sense of 
these parameter choices is to ask how many susceptible 
individuals must be in a cell before a single zombie has 
a higher rate for biting a human than transitioning to a 
neighboring cell. For our choice of parameters, this gives 

7V/3 = 4/i => N ~ 102 . (46) 

This corresponds to a low population density of 
~ 11 people/km 2 , again agreeing with our intuition. All 
of our parameter choices are summarized in Table [I] 

C. Simulation Details 

To effectively simulate an outbreak at this scale, we 
employed the Next Reaction Method of [6j. We main¬ 
tained a priority queue of all possible reactions, assign¬ 
ing each the time at which the reaction would take place, 
an exponentially distributed random number with scale 
set by the rate for the reaction. At each time step of 
the simulation, we popped the next reaction off of the 
queue, and updated the state of the relevant squares on 
our grid. Whenever population counts changed, we of 
course needed to update the times for the reactions that 
depend on those population counts. This method re¬ 
mained efficient for simulating the entire US. However, 
at late times a large amount of simulation time was spent 
simulating the diffusion of the zombies back and forth be¬ 
tween highly populated states. We could have achieved 
additional computational efficiency by adopting the time 
dependent propensity function approach of Fu et al. [5]. 

D. Results 

With the simulation in place, we are now in a position 
to simulate a full scale zombie outbreak. We first consider 
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an outbreak that began with one in every million indi¬ 
viduals starting in the Exposed (E) state in the United 
States. For a single instance the overall populations are 
shown in Figure [lO] This looks similar to the analytical 
outbreaks we saw in Figure [lj but with a steeper rate 
of initial infection and some slight perturbations to the 
curves. The total population curves however hide most of 
the interesting features. In Figure [TT] we attempt to give 
a sense of how this outbreak evolves, showing the state 
of the United States at various times after the outbreak 
begins. 



FIG. 10. The S (blue), Z (red), R (black), and E (green) pop¬ 
ulations as a function of time for a full scale zombie outbreak 
in the continental United States starting with one in every 
million people infected (color online). The exposed popula¬ 
tion (E) has been magnified by a factor of 100. 

As you can see, for the parameters we chose, most of 
the United States population has been turned into zom¬ 
bies by the first week, while the geographic map does 
not necessarily seem all that compelling. In the early 
stages of the outbreak, while the population is roughly 
homogeneous, the zombie plague spreads out in roughly 
uniform circles, where the speed of the infection is tied to 
the local population density. Infestations on the coasts, 
with their higher population density, have spread far¬ 
ther than those near the center of the country. After 
several weeks, the map exhibits stronger anisotropy, as 
we spread over larger geographical areas and the zombie 
front is influenced by large inhomogeneities in population 
density. After four weeks, much of the United States has 
fallen, but it takes a very long time for the zombies to 
diffuse and capture the remaining portions of the United 
States. Even four months in, remote areas of Montana 
and Nevada remain zombie free. 

To investigate the geographical characteristics of the 
outbreak, we must move beyond a single instance of an 
outbreak and study how different regions are affected in 
an ensemble of outbreaks. If it takes a month to develop 
and distribute an effective vaccine (or an effective strat¬ 
egy for zombie decapitation), what regions should one lo¬ 
cate the zombie-fighting headquarters? We ran 7,000 dif¬ 
ferent 28-day zombie outbreaks in the continental United 
States starting with a single individual. A single instance 
of one of these outbreaks originating in New York City 



(e) 3 Weeks (f) 4 Weeks 



(g) 2 Months (h) 4 Months 


FIG. 11. Simulation of a zombie outbreak in the continental 
United States. Initially one in every million individuals was 
infected at random. Results are shown above at (a) one day, 
(b) two days, (c) one week, (d) two weeks, (e) three weeks, 
(f) four weeks, and (g) two months after the outbreak begins. 
Shown here are the population of susceptible individuals (S) 
in blue, scaled logarithmically, zombies in red and removed in 
green (color online). All three channels are superimposed. A 
movie version of this outbreak is available in the supplemental 
materials online. 


is shown in Figure [T2| 

By averaging over all of these runs, we can start to 
build a zombie susceptibility map, as shown in Figure [13) 
In the top plot, we show the probability that the given 
cell is overrun by zombies after seven days. Here you 
can clearly see that there are certain regions - those sur¬ 
rounding populous metropolitan areas - that are at a 
greater risk. This is partly because those regions have 
lots of individuals who could potential serve as patient 
zero, and partly due to the rapid spread of zombies in 
those areas. In the bottom plot, we plot the probability 
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FIG. 12. Status of the United States 28 days after an out¬ 
break that started in New York City. Here blue represents 
humans, red represents zombies and green represents dead 
zombies (color online). The three color channels have been 
laid on top of one another. 


that the cell is overrun, but at the 28 day mark. 

After 28 days, it is not the largest metropolitan ar¬ 
eas that suffer the greatest risk, but the regions located 
between large metropolitan areas. For instance, in Cali¬ 
fornia it is the region near Bakersfield in the San Joaquin 
Valley that is at the greatest risk as this area will be over¬ 
run by zombies whether they originate in the San Fran¬ 
cisco area or the Los Angeles / San Diego area. The area 
with the greatest one month zombie risk is north eastern 
Pennsylvania, itself being susceptible to outbreaks origi¬ 
nating in any of the large metropolitan areas on the east 
coast. 


VI. CONCLUSION 

Zombies offer a fun framework for introducing many 
modern concepts from epidemiology and critical phenom¬ 
ena. We have described and analyzed various zombie 
models, from one describing deterministic dynamics in a 
well-mixed system to a full scale US epidemic. We have 


given a closed form analytical solution to the well-mixed 
dynamic differential equation model. We compared the 
stochastic dynamics to a comparable density-dependent 
SIR model. We investigated the critical behavior of the 
single person per site two-dimensional square lattice zom¬ 
bie model and demonstrated it is in the percolation uni¬ 
versality class. We ran full scale simulations of a zombie 
epidemic, incorporating each human in the continental 
United States, and discussed the geographical implica¬ 
tions for survival. 

While this work is predicated on a fictional infesta¬ 
tion, one might ask whether there are any phenomena in 
the real world that behave in a manner similar to our 
modeled zombie outbreaks. As noted, the SZR model 
requires that susceptible hosts directly participate in the 
removal of zombie hosts from the infectious population, 
leading to runaway outbreaks as susceptible hosts are 
depleted. One might imagine a similar phenomenon for 
infectious diseases that require medical intervention to 
be suppressed; as medical personnel themselves become 
infected (as has sadly happened to a considerable degree 
during the recent Ebola outbreak in West Africa), they 
become less able to stem the spread of infection. (Med¬ 
ical personnel, however, represent only a small fraction 
of all susceptible hosts, so a refinement to an SZR -type 
model would be required to account for this.) One might 
also imagine SZR-like dynamics in the spread of ideas 
and opinions: a person spreading a controversial opin¬ 
ion in a population, for example, might be able to sway 
some converts, but is also likely to meet resistance and 
counter-arguments, which act to reduce infectivity and 
perhaps ultimately stop the spread. 

We hope our systematic treatment of an imaginary dis¬ 
ease will provide a useful and inspiring teaser for the ex¬ 
citing fields of statistical mechanics, network science, and 
epidemiology. 
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